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Abstract 

Parametric two-electron reduced-density-matrix (p-2RDM) methods have enjoyed much success 
in recent years; the methods have been shown to exhibit accuracies greater than coupled cluster 
with single and double substitutions (CCSD) for both closed- and open-shell ground-state ener- 
gies, properties, geometric parameters, and harmonic frequencies. The class of methods is herein 
discussed within the context of the coupled electron pair approximation (CEPA), and several 
CEPA-like topological factors are presented for use within the p-2RDM framework. The resulting 
p-2RDM/ra methods can be viewed as a density-based generalization of CEPA/n family that are 
numerically very similar to traditional CEPA methodologies. We cite the important distinction 
that the obtained energies represent stationary points, facilitating the efficient evaluation of prop- 
erties and geometric derivatives. The p-2RDM/re formalism is generalized for an equal treatment of 
exclusion-principle- violating (EPV) diagrams that occur in the occupied and virtual spaces. One of 
these general topological factors is shown to be identical to that proposed by Kollmar [C. Kollmar, 
J. Chem. Phys. 125, 084108 (2006)], derived in an effort to approximately enforce the D, Q, and 
G conditions for iV-representability in his size-extensive density matrix functional. 
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I. INTRODUCTION 



It has long been understood that the ground-state energy for a many-electron system can 
be expressed as a functional of the two-electron reduced-density-matrix (2-RDM) l|, |2| and 
thus the energy may be determined without knowledge of the iV-electron wave function. The 
direct determination of the 2-RDM cannot be accomplished without imposing the so-called 



iV-representability conditions [l[ 3, Q]: constraints that guarantee that the 2-RDM corre- 
sponds to a realistic iV-body density. Two general classes of methods for the determination 
of the 2-RDM without knowledge of the many electron wave function have emerged. The 
2-RDM may be determined (i) variationally, whereby the energy is minimized with respect 
to the elements of the 2-RDM under the constraint that the eigenvalues of the matrix re- 



main non-negative 



4ll|. or (ii) non- variational 



contracted Schrodinger equation (ACSE) 12|-[22|]. In the non- variational formalism, the 



y, via the solution of the anti-Hermitian 



iV-representability conditions take the form of the cumulant reconstruction of the 3- RDM. 

The D, Q, and G conditions for iV-representability have been approximately incorporated 
into a size-extensive functional of the 2-RDM as parametrized by a set of single and dou- 
ble excitation coefficients 



23 



30]. This parametric approach to variational 2-RDM theory 



exhibits accuracies superior to CCSD in terms of energies, properties, geometric parame- 
ters, and harmonic frequencies. The methodology has also been extended to include general 
forms for the equivalent treatment of both open- and closed-shell systems and for use within 
local correlation approximations suitable for the treatment of very large molecules 29J]. Re- 
cently, an alternate parametrization has been developed developed by one of the authors 
that stresses the equal treatment of occupied and virtual spaces exhibitingthe accuracy 
of CCSD with perturbative triple excitations (CCSD(T)) at equilibrium 27J, and unlike 
CCSD(T), the quality of the 2-RDM solution does not degrade when stretching a single 
chemical bond. 

The cou pled electron pair approximation (CEPA) 3lM35| and coupled pair functional 



(CPF) [36 



371 ] enjoyed much success in the 1970s, but fell out of use with the advent of 



efficient vectorized coupled-cluster algorithms. By disregarding some diagrams describing 
disconnected triple and quadruple excitation amplitudes, coupled pair theories can be im- 
plemented at a cost that is slightly less than CCSD. The methods are therefore often viewed 
as at best an approximation to CCSD. Both early and recent numerical examples, however, 
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have demonstrated that CEPA is no less accurate for single-reference problems. Although 
formally less complete than CCSD, coupled pair methods have many nice properties that 



have made them the subject of considerable interest lately [38H41L l44j . |45[. Provided one 
works within a local orbital basis, CEPA is size-extensive. CEPA is conceptually simpler 
than CCSD, and while formally exhibiting the same scaling, may be more suitable for large- 
scale parallelization. It has recently been demonstrated that the CEPA/1 variant yields 
thermochemical properties intermediate in quality between CCSD and CCSD(T) 40|. Fur- 
thermore, the CEPA methods can be expressed within the framework of the CPF and the 



topological matrix originally proposed by Ahlrichs 



implying that the CEPA variants 



may 
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be incorporated into a CI framework via a diagonal shift to the Hamiltonian matrix 



41] . As a result, efficient CI algorithms can be utilized to solve the CEPA equations. 



It is clear that there are connections between the topological matrix of Cl-driven CEPA 
and the topological factor that lends the parametric 2-RDM method its size-extensivity 
properties. In this paper, we derive CEPA-like overlap equations for the parametric 2-RDM 
method to elucidate the connection between CEPA and 2-RDM methods. In this formalism, 
the 2-RDM method appears very similar to CEPA with two key exceptions. First, unlike 
traditional CEPA methods, parametric 2-RDM methods account for the so-called exclusion- 
principle- violating (EPV) diagrams that occur in the virtual space; these considerations are a 
consequence of the balanced treatment of particles and holes that emerges when considering 
the D, Q, and G conditions for iV-representability. Second, the overlap equations for the 
parametric 2-RDM method contain an additional term that renders the solution stationary 
in all of its variables. Accordingly, the determination of density matrices and thus properties 
and the evaluation of geometric derivatives are all greatly simplified. In fact, the derivative 
of the energy with respect to an arbitrary nuclear coordinate requires only the derivatives of 
the 1- and 2-electron integrals, which may be evaluated either analytically or numerically. 

From explicit relations between the non- variational CEPA and the variational parametric 
2-RDM methods based on the observations above, we derive topological factors for the 
parametric 2-RDM methods that correspond to the CEPA/n variants (with n = 1,2,3), 
as well as new factors that correspond to the same type of hierarchy, but with a balanced 
treatment of occupied and virtual spaces. The topological factor proposed by Kollmar 3] 
is shown to be one that corresponds to variant 1 of a CEPA theory that accounts for EPV 
diagrams in the virtual space. We demonstrate the similarities of the two methodologies 
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numerically with applications to bond stretches and geometry optimizations for several small 
molecules. We also demonstrate the necessity of the balanced description of particles and 
holes with a bond stretch for the CH radical. In this difficult case, most single-reference 
theories exhibit a qualitatively unphysical hump in the potential energy surface; accounting 
for EPV diagrams in the virtual space alleviates this problem. Finally, it is well known that 
the CEPA methodologies (with the expectation of the CEPA/0 variant) are not invariant 
to unitary transformations of the occupied orbitals. We investigate the numerical behavior 
of the 2-RDM method with respect to the same types of orbital rotations. The 2-RDM 
method is shown to vary slightly with the choice of orbitals, but this variance is insignificant 
compared to the total correlation energy. The 2-RDM method is shown to be rigorously 
size-extensive for noninteracting two-electron and two-hole systems when these systems are 
described by a basis of localized molecular orbitals. 

The key relations of this paper between the non-variational CEPA and the variational 
parametric 2-RDM methods were first presented by DePrince in his Ph.D. thesis at The 
University of Chicago in 2009 42] . Similar results have a ppe ared in a recent paper, published 
to the web on August 18, 2010, by Neese and Kollmar 43[. To facilitate circulation of our 
work, we publish the present paper to the Archives while we finish a more complete version 
of the paper for publication elsewhere. 



II. THEORY 

In Section Hi Al we briefly review parametric 2-RDM theory. Section Hi Bl outlines the cou- 
pled electron pair approximation, and Section Hi CI provides a link between the two method- 
ologies. 



A. Parametric 2-RDM methods 



To most easily elucidate the connections between parametric 2-RDM and CEPA theories, 
we limit our discussion at this point to the configuration interaction wave function with only 
double excitations: 

|*) = Co|* ) + £ c ?in 6 )> (1) 

a<b 

i<j 
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where |\l/o) represents the reference wave function, represents a doubly substituted 

configuration in which occupied spin-orbitals i and j have been replaced by virtual spin- 
orbitals a and b, and the set of coefficients {c , cfj*} represents the respective CI expansion 
coefficients for these configurations. For the normalized wave function, the well-known lack of 
size-extensivity associated with truncated CI methods is attributed to those terms in which 
the excited determinants interact with the reference configuration. Size-extensivity may be 
restored by the introduction of a purely connected generalized normalization coefficient into 
the corresponding energy expression, 

E c = Y,^o\Hm>)c%c £ + £ - EolV*)^, (2) 

a<b a<b c<d 

i<j i<j k<l 

where the generalized leading coefficient, c |f, is defined according to Kollmar's definition 



23) as 



c j? = (i-Ei4?r v$f) i/2 - (3) 

c<d 
k<l 

The 8-index topological matrix, 2 /^ d , interpolates between iV-representable the (but not 
size-extensive) CID solution, recovered by setting all 2 f^t equal to unity, and the size- 
extensive (but not iV-representable) CEPA/0 solution, obtained by setting all Vij^f equal 
to zero. It is true that the CEPA/0 parametrization will restore size-extensivity, but we 
arrive at a more intelligent choice by recognizing that the iV-representability of the as- 
sociated 2-RDM is strictly dependent upon the form of the topological factor. We will 
choose this factor such that it will enforce, at least approximately, the known two-particle 
iV-representability conditions, the so-called D, Q, and G conditions. With these considera- 



tions, Kollmar proposed in Ref. 23[ the following factor: 



cabcd rpkl , rpcd rpkl rpcd I a\ 

Jijkl - *ij + *ab - *ij *ab > l 4 J 



1 



F;l=^(S ps + S pt + 5 qs + 5 qt ), (5) 

where the delta function 5 pq is zero when the spatial component of orbitals p and q are 
disjoint and one otherwise. By replacing c ^ fc in Eq. ([2]) with its definition in Eq. ([3]), the 
correlation energy may be determined via an unconstrained minimization. 
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B. The coupled electron pair approximation 



Beginning with the CID energy functional, a set of non-linear equations may be obtained 



by enforcing the stationary condition, dE c jdcf- = 0, to obtain 



„ab 

"ij 

,,u6 



Ec = E— <*?l#l*o>, (6) 

a<b °0 

d*E e = Co(^\H\%) + Y,(W$\H ~ EomlKt (7) 

c<d 
k<l 

Making a change of variables, = c^/c , yields the CID overlap equations in intermediate 
normalization: 

E c = Y,b%(^\mo), (8) 

a<b 

i<j 

E c b% = (*#|ff|¥ ) + E<*?l# ~ EoinM. (9) 

c<d 
k<l 

The left hand side of Eq. is clearly quadratic; the size-extensivity problem in this 
representation amounts to the lack of a complementary quadratic term on the right hand 
side of this equation. Such a term describes the interaction between all doubly and quadruply 
excited configurations: 

= (*£|£|¥ > + J2(^\H -E - E e \9i)US + (*#|&|*g), (10) 

c<d 
k<l 

iabcd\ 



termediately normalized CI expansion coefficients, b^f. The coupled-cluster with doubles 
(CCD) equations are recovered by approximating the coefficients of the quadruples, b^j: d 



where \^q) includes all quadruply excited configurations, (^"j^f), and their respective in- 

'abcd 
Hjkl 

Ti^inlri /ii nnr n i~\ ■ r Kin /Nil nnvnrxlnn I, 

ijkl 

as an antisymmetric sum of products of doubles coefficients as suggested by second order 
perturbation theory. The CEPA methods imply a simpler relationship between the double 
and quadruple coefficients by taking only the leading term of this sum: 

U$£ = 6?6g- (11) 

Inserting Eq. ( TTTj) into Eq. ( TlOT) yields 

= «|£|*o> + -Eo- S e |*g>6g + Y,(^\Hm-H)bpl1, (12) 

c<d c<d 
k<l k<l 



and the last term may be reexpressed equivalently according to Slater's rules to give 

= (9/$\H\* ) + £<*#|# -Eo- E e \*$Vg + £<*o|#|*g>$&g, (13) 

c<d c<d 
k<l k<l 

or 

= «|£|*o> + J2(^\H -Eo- £ e |*g>fig + EJ%. (14) 

c<d 
k<l 

We have arrived at the simplest CEPA approximation, denoted CEPA/0: 

= (*£|#|¥o> + £<*#|# " 3>l*w>&£i- (15) 

The CEPA/0 approximation is a naive one in that we have unintentionally included the 
effects of a large number of unphysical terms. Equation (fT5l) deteriorates whenever the 
indices {ij} and {kl} (or {a&} and {cd}) have any coincidences. Such instances imply 
multiple excitations out of (or into) the same orbitals twice and are therefore referred to 
as exclusion-principle-violating (EPV) terms. The remainder of the CEPA approximations 
differ only in how they account for EPV terms. 

By defining a diagonal shift, A?-, we may write general equations for all of the CEPA 
variants that improve upon CEPA/0, 

o = (*$\h\* ) + E<*?l# ~ E »- A#|*g>&& (16) 

c<d 
k<l 

where we can easily see that A^ b is equal to — E c for CID and zero for CEPA/0. The simplest 
improvement upon CEPA/0, called CEPA/2, removes only those EPV terms of the form 
bfjbij - ; the diagonal shift for CEPA/2 is defined as the pair energy, e^, 

^ = ^ = T.^m-m^o). (it) 

c<d 

Table H] lists the definitions of the diagonal shifts for the CID and CEPA/n (n = 0, 1,2,3) 
methods. CEPA/3 removes all EPV diagrams in the occupied space, and the CEPA/1 shift 
can be viewed as the average of those corresponding to CEPA/2 and CEPA/3. In general, 
the number of EPV diagrams removed by each flavor of CEPA is CEPA/3 > CEPA/1 > 
CEPA/2 > CEPA/0, and as such the correlation energy is lowest for CEPA/0 and highest 
for CEPA/3. We may express in a more suggestive form by incorporating idea of the 



7 



TABLE I: Diagonal shifts that define the CID and CEPA equations. The various CEPA shifts lend 
size-extensivity to CID while removing varying degrees of unphysical exclusion principle violating 
(EPV) terms from the overlap equations. 



Method A# 

CID -E c 

CEPA/0 

CEPA/1 ^ k (e lk + e jk ) 

CEPA/2 eij 

CEPA/3 £ fe (e ife + e jk ) - e* 



topological factor into Eq. ( TTTj) . and symmetrizing the factor with respect to the exchange 
of orbitals i and j or k and I, 

^ = T.inimo^ 2 f$ k f, (is) 

c<d 
k<l 

where 

2 ftjki = -^{Sikbji + 5 jk Su). (19) 
It is immediately clear that one could define a topological factor corresponding for each 
of the CEPA/n variations; these factors, symmetrized with respect to orbital exchange, are 
presented in Table ILT1 We will show in the next section that these factors can be incorporated 
into the parametric 2-RDM energy functional to yield a family of density-based CEPA-like 
methods, which we call p-2RDM/n. 

C. Density-based CEPA 

In this section, we illustrate the connection between the parametric 2-RDM method and 
the CEPA/n family of equations. By enforcing the stationary condition on Eq. (j2J), we 
obtain the following system of coupled nonlinear equations that define the parametric 2- 
RDM energy and excitation coefficients: 

„ab 

^ = E-4f(^o|^|^), (20) 

a<b C 0,ij 
i<j 

8 



TABLE II: Symmetrized topological factors corresponding to the CEPA/ra family of methods. Note 
that EPV diagrams are neglected in virtual space. When incorporated into the parametric 2-RDM 
method, the methods are denoted p-2RDM/n. 



Method /$f 

CID 1 

CEPA/0 

CEPA/1 i(5. k + s jl + 6 u + 5 jk ) 

CEPA/2 + 

CEPA/3 \{$ik + Sji + Sii + 5jk - (SikSji + 5ii5jk)) 



o = (*g|ff|* >cSf„ - 2 /g? (*o|^|^>-4f + - 3>l*?><#- (2i) 

a<6 C 0,ij a<fe 

We can make the transformation bfj = c?- / c to obtain a set of equations that is identical 
to the CEPA family of equations given in Eq. ( |T6l) with a diagonal shift that is defined as 

/ c cd \ icd 

A- 6 = E & h 2 fm «\h\*o) + E - 1 7S(n 6 l^l^)- ( 22 ) 

c<d c<d \ C 0,ij / %' 

k<l k<l 

Again, from the perspective of coupled pair theories, A^ can be interpreted as an approxi- 
mation of the effects of higher excitations neglected in the CI expansion. The second term 
in Eq. fl22|) can be seen as one that renders the energy and true minimum and thus that 
the corresponding density matrix is a stationary solution to Eq. fl5]). Assuming that this 
term is sufficiently small and may be ignored, the connection between parametric 2-RDM 
methods and the CEPA approximations is obvious. This assumption is not unreasonable 
for well-behaved systems; the term is exactly zero in the both the CID and CEPA/0 limits. 
By replacing the topological factor in Eq. ([2]) or fl22|) with those denned in Table [Til we 
obtain a family of methods with stationary solutions that yield results that are numerically 
very similar to traditional CEPA/n implementations. We term this family of density-based 
CEPA-like methods parametric 2-RDM/n or p-2RDM/n methods. 

Returning to the original formulation of the parametric 2-RDM method, the Kollmar 
topological factor, which we will call K, can be viewed as one that is very similar to CEPA/1 
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with the important distinction that it accounts for EPV diagrams in the virtual space. This 
balanced description of occupied and virtual spaces is central to reduced-density-matrix 
theory (consider the complementary D, Q, and G conditions for iV-representability). We 
may define a family of density-based CEPA-like topological factors to be incorporated into 
the parametric 2-RDM formalism by (i) accounting for EPV diagrams in the virtual space 
and (ii) symmetrizing each factor with respect to the exchange of orbitals % and j, k and I, 
a and b, or c and d. These factors are defined in Table IHIl These new topological matrices 
yield a family of improved density-based CEPA methods and are labeled p-2RDM'/n (for 

TABLE III: Symmetrized topological factors with a balanced description of the occupied and 
virtual spaces. Each factor, with the exception of CEPA/0, yields the exact correlation energy for 
two-electron systems. Each topological factor is defined as a combination of tensors corresponding 
to the occupied and virtual spaces: f^f = F-f + F^ — F-fF?£. 



Method F* 

CID 1 

CEPA/0 

p-2RDM'/l l(s ps + 6pt + 6gs + 6qt ) 

p-2RDM'/2 \ (SpsSgt + 5 pt 5 qs ) 

p-2RDM'/3 %(5p3+5pt + 

8qs + Sgt) 



n = 1,2,3). Each p-2RDM'/n variant may be implemented at a cost that is comparable 
to other two-electron theories. It should be noted that each factor, with the exception 
of CEPA/0, recovers the exact correlation energy in the two-particle (or two- hole) limit. 
Furthermore, when incorporated in the parametric 2-RDM formalism, which is Hermitian, 
geometry optimizations, the determination of density matrices, and the evaluation of one- 
and two-electron properties are all greatly simplified. 

To this point we have not considered the effects of single excitations on size-extensivity 
and EPV diagrams. No clear consensus for the treatment of single excitations in coupled-pair 
formalisms is present in the literature. For this reason, we choose to treat single excitations 
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in each of the p-2RDM/n and p-2RDM'/n methods as is described in Ref. 28|. We have 
the generalized normalization coefficient, c Q , defined as 



„ u.u _ (■, ST^ \ n c\2 1 rabcc | cd 1 2 2 f abcd\ 1/2 /r>o\ 

C n ^- — ^-L 2^\ C k\ Jijkk 2^\ C kl\ Jijkl) J l zo J 



aft 

c c<d 
k k<l 



where 1 ff^ is either defined as 



VgS = l-(l-**)(l-*i*)> ( 24 ) 
in the p-2RDM/n formalism, or 

VgS = 1 - (1 - <W(1 - S jk )(l - 5 ac )(l - <U (25) 

in the improved p-2RDM'/n formalism. The value of as given by Eq. (12411 is unity 

unless the spatial component of the occupied orbitals i and j are disjoint with k; this is 
the treatment chosen to coincide the best with existing CEPA theories. The value of 1 
as given by Eq. (1241) is unity unless the spatial components of the occupied orbitals i 
and j are disjoint with k and the spatial components of the virtual orbitals a and b are 
disjoint with c. From the perspective of coupled pair theories, our singles topological factors 
approximate the inclusion of disconnected triple excitations while removing EPV diagrams 
in the occupied space or occupied and virtual spaces. We note that this treatment of single 
excitations removes all EPV diagrams in either the occupied space or occupied and virtual 
spaces arising from single excitations and is thus most similar to the CEPA/3 treatment of 
single excitations. 



III. DISCUSSION 



CEPA/n calculations were performed using the Molpro electronic structure package 46]. 



The closed-shell p-2RDM/n calculations were performed using our code implemented within 
the PSI3 ah initio electronic structure package [30I . I^rj ] . All open-shell parametric 2-RDM 
calculations were performed with a separate code with all 1- and 2-electron integrals obtained 
from the GAMESS electronic structure package j^. The topological factors for each p- 
2RDM/n and p-2RDM'/n variant are presented in Tables HT1 and [TTTl respectively; note 
that, as in the Molpro implementation of CEPA/ri, EPV diagrams in the virtual space are 
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neglected for p-2RDM/n calculations. Single excitations are treated in the p-2RDM/n and 
p-2RDM'/n methods as described in Ref. 28 1. 

Figure [T] illustrates the CEPA/n and p-2RDM/n potential energy surfaces for a single O- 
H bon stretch for H 2 in a cc-pVDZ basis set. Clearly the p-2RDM/n and CEPA/n methods 



FIG. 1: Potential energy curve for a single O-H bond stretch for H2O in a cc-pVDZ basis set with 
one core orbital frozen. Curves are presented for the CEPA and p-2RDM variants 1, 2, and 3. 
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perform identically in all regions of the potential energy curve. The largest deviation between 
the two families is 0.95 milli-Hartrees (mH), occurring for p-2RDM/2 at a bond length of 
2.2 A. Molpro CEPA/2 calculations did not converge beyond 2.2 A. The largest discrepancy 
with CEPA/3 is only 0.47 mH, occurring at 3.0 A. Figure [2] illustrates similar trends for the 
NH 2 -H bond stretch. CEPA/n and p-2RDM/n results are indistinguishable over the range 
of reported bond lengths. 

The p-2RDM/n and CEPA/n formalisms were also applied geometry optimizations and 
harmonic frequency analysis for H2O and CO2 in a cc-pVDZ basis set. The optimized bond 
lengths and energies for CO2 as computed by the CEPA/n and p-2RDM/n methods are 
listed in Table IIV( in general we can see that the optimal bond length contracts with a 
more rigorous treatment of EPV diagrams. The CEPA/3 and p-2RDM/3 results in Table 
IIVI are indistinguishable. The discrepancies that arise between the variant 2 and 1 numbers 
are due to the difference in the treatment of single excitations in CEPA/n and p-2RDM/n 
theories. The p-2RDM/n methods remove all EPV diagrams in the occupied space due 
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FIG. 2: Potential energy curve for a single N-H bond stretch for NH3 in a cc-pVDZ basis set with 
one core orbital frozen. The bond length for one hydrogen is increased while holding all other 
bonds and angles constant. Curves are presented for the CEPA and p-2RDM variants 1, 2, and 3. 
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TABLE IV: Optimized geometries and energies for CO2 in a cc-pVDZ basis set. Energies and bond 
lengths are given in Hartrees and A, respectively. CEPA/n and p-2RDM/n methods yield nearly 
identical results for both the bond length and the minimum energy. Core orbitals are restricted to 
be occupied. 



Variant 


CEPA 




P-2RDM 




Energy 


rco 


Energy 


rco 


1 


-188.1347 


1.1713 


-188.1341 


1.1707 


2 


-188.1440 


1.1743 


-188.1424 


1.1729 


3 


-188.1263 


1.1688 


-188.1263 


1.1688 



to single excitations whereas CEPA/3 removes more than CEPA/1, which removes more 
than CEPA/2. Accordingly, we observe the greatest disparities between the variant 2 val- 
ues. Tables |V] and I VI I present the optimized energies, geometric parameters, and harmonic 
frequencies as computed by the CEPA/n and p-2RDM/n methods for the H2O molecule. 
Geometric parameters are identical within each variant, with a difference in bond angle of 



13 



TABLE V: Optimized geometries and energies for H 2 in a cc-pVDZ basis set. Energies, bond 
lengths, and angles are given in Hartrees, A, and degrees, respectively. CEPA/n and p-2RDM/n 
methods yield identical results for the geometric parameters and the minimum energy. Core orbitals 
are restricted to be occupied. The experimentally obtained values for roH and aHOH are 0.9578 
A and 104.4776 degrees, respectively. 







CEPA 






P-2RDM 




Variant 


Energy 


roH 


a HOH 


Energy 


roH 


a HOH 


1 


-76.2382 


0.9652 


102.11 


-76.2382 


0.9652 


102.12 


2 


-76.2406 


0.9663 


101.99 


-76.2406 


0.9663 


102.00 


3 


-76.2360 


0.9642 


102.22 


-76.2360 


0.9642 


102.22 


TABLE VI: Harmonic frequencies in wavenumbers, cm" 


1 , for H2O in a cc-pVDZ basis set. 


CEPA/n 


and p-2RMD/n methods yield are nearly indistinguishable, with the largest discrepancies being 


1.9 cm 1 


for the symmetric 


and asymmetric 


stretches 


as described by variant 2. 








CEPA 






P-2RDM 




Variant 


ai 


ai 


b 2 


ai 


ai 


b 2 


1 


3837.1 


1695.3 


3937.8 


3838.2 


1696.4 


3938.7 


2 


3815.9 


1691.8 


3917.8 


3817.8 


1692.9 


3919.7 


3 


3855.9 


1698.5 


3955.3 


3855.5 


1699.4 


3955.5 



only 0.01 degrees for variants 1 and 2. As was observed with CO2, bond lengths slightly 
contract with a more rigorous treatment of EPV diagrams. Harmonic frequencies are nearly 
identical between the two methods, with the largest difference between CEPA and p-2RDM 
being less than two wavenumbers. That larger differences exist between geometric parame- 
ters determined by CEPA and p-2RDM methods in the case of CO2 rather than H 2 is not 
a surprise. The treatment of single excitations is directly connected to the description of 
disconnected triple excitations, and the emergence of discrepancies between the two methods 
is most likely due to the growing importance of triple excitations for C0 2 as compared to 
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H 2 0. 

We next investigate the importance of the equal treatment of EPV diagrams arising in 
the occupied and virtual spaces. In CEPA methodologies, the EPV terms for the virtual 
space are generally ignored because they are far fewer in number than those occurring in the 
occupied space. This omission is simply intended to increase computational efficiency There 
are certain situations, however, where the neglect of virtual space EPV terms may cause the 
CEPA methods to qualitatively fail. Figure |3] illustrates such a point. We calculated the 

FIG. 3: Potential energy curve for the bond stretch for the CH radical in a cc-pVDZ basis set with 
one core orbital frozen. CCSD and p-2RDM/3 both exhibit and unphysical hump at long bond 
lengths. p-2RDM'/3 has no such feature. 
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potential energy curve for the CH radical with CCSD and the p-2RDM/3 and p-2RDM'/3 
methods using the topological factors given in Tables HT1 and 1111} respectively. At around 2.8 
A an unphysical hump occurs in the CCSD curve. p-2RDM/3 neglects the virtual space EPV 
diagrams and as a result develops a similar hump around 3.1 A. The p-2RDM'/3 method 
never displays any unphysical behavior. 

It is well known that the pair energies associated with the EPV diagrams are not invariant 



to rotations among the occupied orbitals 



49] . This variance is a consequence of the partial 



nature of the summations that define pair energies; the indices involved to not span the 
entirety of the Hilbert space, and the pair contribution to the energy is thus not invariant 
to unitary transformations. Accordingly, size-extensivity can only rigorously be achieved in 
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these coupled-pair theories with the use of a localized molecular orbital basis. The close 
relationship of the parametric 2-RDM method to these theories suggests that similar defi- 
ciencies may exist within the present formulation of the 2-RDM method. The numerical 
size-extensivity of the parametric 2-RDM method was demonstrated by the authors for a 
series of infinitely separated He atoms in Ref. 24J. We revisit this system, illustrating in Fig. 
Ill II the size-extensivity and orbital invariance properties (or lack thereof) for the paramet- 
ric 2-RDM method. We choose the Kollmar parametrization, K (also denoted p-2RDM'/l 
presently), as the representative example. We treat an increasing number of helium atoms 



FIG. 4: Energy per helium atom as a function of the number of He atoms for a system of non- 
interacting helium atoms. The 2-RDM method is numerically size-extensive utilizing both canonical 
and localized molecular orbitals. Size-extensivity is absolutely rigorous in the local orbital basis. 
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at effectively infinite separation to illustrate the size-extensivity of the 2-RDM method in 
a basis of canonical and localized orbitals. The He atoms are situated on a line with an 
interatomic distance of 200 A; the atomic orbitals are represented by an Ahlrichs double-zeta 



basis set 



50j. For the localized-orbital calculations, occupied and virtual orbitals were lo 



calized separately according to the Boys localization criterion [51| . The energy per He atom 
for a size-extensive method should not vary with the number of He atoms. We see that 
the 2-RDM does display this characteristic when utilizing either localized and delocalized 
canonical orbitals; regardless of the "bumps" in the canonical basis, the 2-RDM results do 
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not display any systematic dependence upon system size. The energies obtained in the local 
basis, however, are rigorously size-extensive, with no numerical deviation in the energy per 
He atom for all system sizes. We have also presented the energies obtained from the David- 
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3, 



son correction [52J, denoted CISD+Q, and the renormalized Davidson correction 
denoted CISD+R(Q). Clearly, neither the Davidson corrected nor renormalized Davidson 
corrected energies are rigorously size-extensive, with both results exhibiting a clear and 
systematic dependence on system size. 

These results unfortunately demonstrate that the 2-RDM energy is not strictly invariant 
to unitary transformations among the occupied orbitals, as is the case with traditional 
CEPA methodologies. Perhaps even more unfortunately, this dependence also extends to 
the virtual space. The dependence upon the choice of the virtual orbitals is a consequence 
of the symmetry properties of the topological factor in the occupied and virtual spaces; 
EPV diagrams are removed from not only the occupied space, as is the case in CEPA, but 
from the virtual space as well. Figure |5] illustrates the deviation of the energy obtained 
from calculations with localized orbitals from those performed using canonical Hartree-Fock 
orbitals for the H-F bond stretch in (a) cc-pVDZ and (b) cc-pVTZ basis sets. We present 
results for two cases within each basis set: (i) calculations in which only the occupied orbitals 
are localized and (ii) calculations in which we localize both the occupied and virtual orbitals 
separately. Clearly, the 2-RDM method is not invariant to unitary transformations in either 
the occupied or virtual subspaces. Importantly, the method's variance with respect to the 
occupied space is nearly constant at all bond lengths and across basis sets. We note that the 
sign of the percent change for local occupied orbitals changes between basis sets, meaning 
that the energy with localized orbitals was lower than the canonical case in the cc-pVDZ 
basis but higher in the cc-pVTZ basis. While the relative magnitudes of the changes are very 
similar, the change in sign implies that we cannot assume a systematic change in energy 
for arbitrary systems and basis sets when utilizing localized versus delocalized canonical 
orbitals. The variance with respect to the virtual space, while being fairly uniform across 
basis sets, is highly dependent upon the geometry of the system. Local virtual orbitals are 
in general very difficult to determine; minima are often local in nature in the localization 
function, and the orbitals themselves may not vary smoothly with nuclear coordinates. It is 
unclear whether the strong dependence of the energy on the virtual space is a consequence of 
these difficulties or some other systematic deficiency inherent to the p-2RDM'/n methods. 
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FIG. 5: The percentage change in the correlation energy, 
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100% when utilizing 



localized molecular orbitals rather than canonical Hartree-Fock orbitals for the H-F bond stretch 
in (a) cc-pVDZ and (b) cc-pVTZ basis sets. 
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We further investigate the dependence of the 2-RDM energy with respect to variations 
in the virtual space with a size-extensivity example similar to the infinitely separated two- 
electron systems treated above; we here investigate the size-extensivity of the 2-RDM method 
of infinitely separated two- hole systems. As in the He example above, a size-extensive 
double-excitation theory should yield the exact result regardless of system size. Figure Hill 
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illustrates the correlation energy for a system of non-interacting HF molecules in a minimal 

FIG. 6: Energy per HF molecule as a function of the number of HF molecule for a system of 
non- interacting HF molecules. The 2- RDM method is only rigorously size-extensive when using a 
basis of localize occupied and virtual molecular orbitals. 




Number of HF molecules 



ST0-6G basis set. The molecules lie parallel to one another in a line with a distance 
between each center of mass of 1000 A. The H-F bond length is taken as the experimentally 
determined value given in the computational chemistry comparison and benchmark database 



(CCCBD) 55]. The 2- RDM method is utilized with either canonical or separately localized 
occupied and virtual orbitals. As expected, CISD is not size-extensive. The 2-RDM method 
yields effectively size-extensive results for both choices of orbitals, but the only exactly size- 
extensive choice is that in which both the local and virtual orbital spaces are localized. Both 
the Davidson and renormalized Davidson-corrected energies exhibit a strong dependence on 
system size and are thus not rigorously size-extensive. 

IV. CONCLUSIONS 

Parametric variational 2-RDM methods are an accurate and efficient alternative to tra- 
ditional ab initio methods. They formally scale the same as CI with single and double 
excitations and may be implemented at a cost that is slightly less than coupled cluster 
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with single and double substitutions. Provided calculations are performed in a local or- 
bital basis, the obtained energies are rigorously size-extensive. The methods have been 
previously been generalized for geometry optimizations, harmonic frequency analysis, the 
treatment of open-shell systems, and local correlation approximations with much success. 
Novel parametrizations, unique from that originally proposed by Kollmar or discussed in 
this paper have been presented that result in accuracies similar CCSD(T). 

The parametric 2-RDM approach represents a fairly young class of methods, and as such 
its relationship to other methods has remained to this point unexplored. For this reason, 
we have drawn connections between parametric 2-RDM methodologies and existing coupled 
electron pair approximation (CEPA) theories. We have derived a set of topological factors 
that correspond to the CEPA/n (n = 1,2,3) family, and the resulting class of methods is 
a density-based generalization of CEPA/n called p-2RDM/n. Extensive numerical studies 
of equilibrium energies, geometries, and harmonic frequencies have shown that the p-2RDM 
methods perform very similarly to their CEPA analogues for a variety of closed-shell systems. 
New topological factors have been derived specifically to account for the exclusion-principle- 
violating (EPV) terms that arise in the virtual space that are ignored by standard CEPA 



methodologies. Malrieu and coworkers 



4l| understood the importance of the balance be- 



tween occupied and virtual EPV diagrams in their self-consistent, size-consistent truncated 
CI ( (SC) 2 -CISD ); their proposed method is most similar in spirit to the p-2RDM'/3 variant 
discussed herein. Another factor, denoted p-2RDM'/l, is in fact identical to that proposed 
by Kollmar [23J]. The proper treatment of the virtual space EPV diagrams is necessary in 
some situations to obtain physically meaningful results, as was demonstrated for the poten- 
tial energy surface for the CH radical. Consideration of the virtual space EPV diagrams 
for the p-2RDM/3 method is necessary to avoid an unphysical hump in the dissociation 
curve. Aside from these numerical arguments, properly treating virtual space EPV dia- 
grams is absolutely essential from the standpoint of density matrix theory in that they are 
required for a balanced treatment of particles and holes. The parametric 2-RDM framework 
for p-2RDM is also quite convenient as compared to the standard overlap equation formu- 
lation of traditional CEPA methodologies; the energies obtained from any of the p-2RDM 
or p-2RDM' variants presented herein are stationary points, facilitating the evaluation of 
geometric derivatives and the direct computation of density matrices and their associated 
one- and two-electron properties. 
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We have discussed in detail the orbital invariance properties of the CEPA and p-2RDM 
methods citing various numerical examples. The 2-RDM method is indeed exactly size- 
extensive, but this claim holds true only under the condition that one utilizes a basis of 
localized molecular orbitals. As such, the 2-RDM method (and thus the p-2RDM and 
p-2RMD' variants) are not rigorously invariant to unitary transformations within orbital 
subspaces. The p-2RDM methods display a dependence on the choice of orbitals for the 
occupied space while the p-2RMD' methods also vary with the choice of the virtual orbitals. 
One may circumvent any ambiguities with respect to the definition of the orbital space by 
always utilizing a basis of local orbitals. Determining these orbitals in the occupied subspace 
is trivial by the Boys localization criterion 5l|, but may prove problematic for the virtual 
space where local orbitals may not be unique and do not necessarily vary smoothly with 
nuclear coordinates. Fortunately, the variance in the correlation energy with the occupied 
orbitals is only a fraction of a percent and at worst on the order of one percent for the virtual 
orbitals. These discrepancies are very small when compared to the percentage of correlation 
energy that is not recovered for any of the standard ab initio methods with respect to the 
exact full CI results. 
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